

###########################################################
##### Haiti elite network project  		          			#####
##### price regressions	                           		#####
##### 2021 mar 03                   									#####
###########################################################


pri3 <- read.dta('01_Data/02_Clean/price_final.dta')
dat2 <- read.csv('01_Data/02_Clean/fam_con_prod_final.csv')

## set number of sims to do for wild bootstrap SEs
boot = 500
boot = 50


#####
## TABLE A4: Summary statistics, product data
#####

dat2$immig_noneur[dat2$immig==0] <- 0

temp <- unique(subset(dat2, select = c('product', 'fam', 'con_final', 
                                       'coup', 'bonw_02_wnind_st', 'syrian', 'pct_0911_dispo', 
                                       'pct_0911_dispo_old')))
temp <- data.table(subset(temp, is.na(temp$product)==F))
setnames(temp, c('pct_0911_dispo', 'pct_0911_dispo_old'), c('share_fam', 'share_firm'))
nown <- temp[,lapply(.SD, FUN = function(x) length(unique(x[is.na(x)==F]))), by = c('product', 'con_final'), .SDcols = 'fam']
setnames(nown, 'fam', 'nown')
temp <- merge(temp, nown, by = c('product', 'con_final'))

temp1 <- temp[,lapply(.SD, FUN = function(x) length(unique(na.omit(x)))), by = 'product', 
              .SDcols = c('con_final', 'fam')]
temp2 <- temp[,lapply(.SD, FUN = function(x) mean(x, na.rm=T)), by = 'product', 
              .SDcols = c('coup', 'bonw_02_wnind_st', 'syrian')]
temp2 <- unique(subset(pri3, select = c('product', 'coup', 'bonw_02_wnind_st', 'immig_noneur',
                                        'cshare')))
temp$share_fam <- ifelse(temp$share_fam=='Inf', NA, temp$share_fam)
temp3 <- temp[,lapply(.SD, FUN = function(x) max(x, na.rm=T)), by = 'product', 
              .SDcols = c('share_firm', 'share_fam')]

tab <- merge(temp1, temp2, by = 'product')

tab <- data.frame(tab)
tab$product <- c('Cigarettes', 'Cola', 'Edible oil', 'Kerosene', 'Evaporated milk',
                     'Corn meal', 'Medicine', 'Furniture', 'Bread', 'Dry peas', 'Fresh fish',
                     'Chicken', 'Rice', 'Sandals', 'Laundry soap', 'Beauty care', 'Raw sugar',
                     'Fabric')
tab <- tab[order(tab$product),]
print(xtable(tab), include.rownames = F)




#####
## TABLE 2: Prices of goods imported by coup participators during autocratic periods
#####

# set rho1-4 to lagged haitian retail prices
pri3$rho1 <- pri3$price_log_lag1; pri3$rho2 <- pri3$price_log_lag2; pri3$rho3 <- pri3$price_log_lag3; pri3$rho4 <- pri3$price_log_lag4

pri3$centXautoc <- pri3$bonw_02_wnind_st*pri3$autoc
pri3$centXquake1 <- pri3$bonw_02_wnind_st*pri3$quake1
pri3$centXautocXcoup <- pri3$bonw_02_wnind_st*pri3$autoc*pri3$coup
pri3$centXquake1Xcoup <- pri3$bonw_02_wnind_st*pri3$quake1*pri3$coup
pri3$coupXhaitieliteconflict <- pri3$coup*pri3$haitieliteconflict_log
pri3$materialconflict_log <- log(pri3$materialconflict + 1)
  


## Model 1
mod1 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             product + date3 - 1, 
           data = pri3)
# clustered SEs
mod1$se <- coeftest(mod1, vcov = cluster.vcov(mod1, pri3$product))[,2]
mod1$vcov1 <- cluster.vcov(mod1, pri3$product)
# wild bootstrap clustered SEs
mod1$vcov2 <- cluster.boot(mod1, pri3$product, R = boot, boot_type = "wild")
mod1$se2 <- coeftest(mod1, vcov = mod1$vcov2)['coupXautoc', 2]
mod1$p2 <- coeftest(mod1, vcov = mod1$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod1$model)])
mean(pri3$coup[inc.obs==T])
# number of clusters
length(unique(pri3$product[inc.obs==T]))
# test for autocorr in error
mod1$bg <- bgtest(mod1)$p.value
mod1$bg_z <- coeftest(bgtest(mod1))['lag(resid)_1',3]


## Model 2
mod2 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 +  
           product + date3 - 1,
           data = pri3)
# clustered SEs
mod2$se <- coeftest(mod2, vcov = cluster.vcov(mod2, pri3$product))[,2]
mod2$vcov1 <- cluster.vcov(mod2, pri3$product)
# wild boot strap clustered SEs
mod2$vcov2 <- cluster.boot(mod2, pri3$product, R = boot, boot_type = "wild")
mod2$se2 <- coeftest(mod2, mod2$vcov2)['coupXautoc', 2]
mod2$p2 <- coeftest(mod2, mod2$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod2$model)])
mean(pri3$coup[inc.obs==T])
# number of clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod2$jsl <- lht(mod2, 'rho1 = 0', vcov. = mod2$vcov)['Pr(>F)'][2,1]

# test for unit root
mod2$coefficients['rho1'] 
mod2$urt <- lht(mod2, 'rho1 = 1', vcov. = mod2$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod2$bg <- bgtest(mod2)$p.value
mod2$bg_z <- coeftest(bgtest(mod2))['lag(resid)_1',3]


## Model 3
mod3 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             product + date3 - 1,
           data = pri3)
# clustered SEs
mod3$se <- coeftest(mod3, vcov = cluster.vcov(mod3, pri3$product))[,2]
mod3$vcov1 <- cluster.vcov(mod3, pri3$product)
# wild bootstrap clustered SEs
mod3$vcov2 <- cluster.boot(mod3, pri3$product, R = boot, boot_type = "wild")
mod3$se2 <- coeftest(mod3, mod3$vcov2)['coupXautoc', 2]
mod3$p2 <- coeftest(mod3, mod3$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod3$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod3$jsl <- lht(mod3, 'rho1 + rho2 + rho3 + rho4 = 0', vcov. = mod3$vcov)['Pr(>F)'][2,1]
# test for unit root
mod3$coefficients['rho1'] + mod3$coefficients['rho2'] + mod3$coefficients['rho3'] + mod3$coefficients['rho4']
mod3$urt <- lht(mod3, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod3$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod3$bg <- bgtest(mod3)$p.value
mod3$bg_z <- coeftest(bgtest(mod3))['lag(resid)_1',3]


## Model 4
mod4 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             product + date3 - 1, 
           data = pri3)
# clustered SEs
mod4$se <- coeftest(mod4, vcov = cluster.vcov(mod4, pri3$product))[,2]
mod4$vcov1 <- cluster.vcov(mod4, pri3$product)
# wild bootstrap clustered SEs
mod4$vcov2 <- cluster.boot(mod4, pri3$product, R = boot, boot_type = "wild")
mod4$se2 <- coeftest(mod4, mod4$vcov2)['coupXautoc', 2]
mod4$p2 <- coeftest(mod4, mod4$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod4$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod4$jsl <- lht(mod4, 'rho1 + rho2 + rho3 + rho4 = 0', vcov. = mod4$vcov)['Pr(>F)'][2,1]
# test for unit root
mod4$coefficients['rho1'] + mod4$coefficients['rho2'] + mod4$coefficients['rho3'] + mod4$coefficients['rho4']
mod4$urt <- lht(mod4, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod4$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod4$bg <- bgtest(mod4)$p.value
mod4$bg_z <- coeftest(bgtest(mod4))['lag(resid)_1',3]


## Model 5
pri3$firmsXautoc <- pri3$num_firms_log*pri3$autoc
pri3$cshareXautoc <- pri3$cshare_log*pri3$autoc
pri3$cigXconf <- ifelse(pri3$product=="Cigarettes", pri3$noneliteconflict_log, 0)
pri3$colaXconf <- ifelse(pri3$product=="Cola", pri3$noneliteconflict_log, 0)
pri3$oilXconf <- ifelse(pri3$product=="Huile comestible", pri3$noneliteconflict_log, 0)
pri3$keroXconf <- ifelse(pri3$product=="Kerosene", pri3$noneliteconflict_log, 0)
pri3$laitXconf <- ifelse(pri3$product=="Lait evapore non sucre", pri3$noneliteconflict_log, 0)
pri3$maisXconf <- ifelse(pri3$product=="Mais moulu", pri3$noneliteconflict_log, 0)
pri3$medXconf <- ifelse(pri3$product=="Medicaments", pri3$noneliteconflict_log, 0)
pri3$meubXconf <- ifelse(pri3$product=="Meubles de salon", pri3$noneliteconflict_log, 0)
pri3$painXconf <- ifelse(pri3$product=="Pain", pri3$noneliteconflict_log, 0)
pri3$poisXconf <- ifelse(pri3$product=="Huile comestible", pri3$noneliteconflict_log, 0)
pri3$poissonXconf <- ifelse(pri3$product=="Poisson frais", pri3$noneliteconflict_log, 0)
pri3$poulXconf <- ifelse(pri3$product=="Poulet", pri3$noneliteconflict_log, 0)
pri3$rizXconf <- ifelse(pri3$product=="Riz", pri3$noneliteconflict_log, 0)
pri3$sandXconf <- ifelse(pri3$product=="Sandales", pri3$noneliteconflict_log, 0)
pri3$savonXconf <- ifelse(pri3$product=="Savon de lessive", pri3$noneliteconflict_log, 0)
pri3$soinsXconf <- ifelse(pri3$product=="Soins esthetiques", pri3$noneliteconflict_log, 0)
pri3$sucreXconf <- ifelse(pri3$product=="Sucre brut", pri3$noneliteconflict_log, 0)
pri3$tisXconf <- ifelse(pri3$product=="Tissus", pri3$noneliteconflict_log, 0)

mod5 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             firmsXautoc + 
             cshareXautoc +
             cigXconf + colaXconf + oilXconf + keroXconf + 
             laitXconf + maisXconf + medXconf + meubXconf +
             painXconf + poissonXconf + poulXconf + 
             rizXconf + sandXconf + savonXconf + soinsXconf + sucreXconf +
           product + date3 - 1, 
           data = pri3)
# clustered SEs
mod5$se <- coeftest(mod5, vcov = cluster.vcov(mod5, pri3$product))[,2]
mod5$vcov1 <- cluster.vcov(mod5, pri3$product)
# wild bootstrap clustered SEs
mod5$vcov2 <- cluster.boot(mod5, pri3$product, R = boot, boot_type = "wild")
mod5$se2 <- coeftest(mod5, mod5$vcov2)['coupXautoc', 2]
mod5$p2 <- coeftest(mod5, mod5$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod5$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod5$jsl <- lht(mod5, 'rho1 + rho2 + rho3 + rho4 = 0', vcov. = mod5$vcov)['Pr(>F)'][2,1]
# test for unit root
mod5$coefficients['rho1'] + mod5$coefficients['rho2'] + mod5$coefficients['rho3'] + mod5$coefficients['rho4']
mod5$urt <- lht(mod5, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod5$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod5$bg <- bgtest(mod5)$p.value


## Model 5a (with centXautoc control)
mod5a <- lm(price_log ~ coupXautoc + coupXquake1 + 
             centXautoc + centXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
              firmsXautoc + 
              cshareXautoc +
              cigXconf + colaXconf + oilXconf + keroXconf + 
              laitXconf + maisXconf + medXconf + meubXconf +
              painXconf + poissonXconf + poulXconf + 
              rizXconf + sandXconf + savonXconf + soinsXconf + sucreXconf +
              product + date3 - 1, 
           data = pri3)
# clustered SEs
mod5a$se <- coeftest(mod5a, vcov = cluster.vcov(mod5a, pri3$product))[,2]
mod5a$vcov1 <- cluster.vcov(mod5a, pri3$product)
# wild bootstrap clustered SEs
mod5a$vcov2 <- cluster.boot(mod5a, pri3$product, R = boot, boot_type = "wild")
mod5a$se2 <- coeftest(mod5a, mod5a$vcov2)['coupXautoc', 2]
mod5a$p2 <- coeftest(mod5a, mod5a$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod5a$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod5a$jsl <- lht(mod5a, 'rho1 + rho2 + rho3 + rho4 = 0', vcov. = mod5a$vcov)['Pr(>F)'][2,1]
# test for unit root
mod5a$coefficients['rho1'] + mod5a$coefficients['rho2'] + mod5a$coefficients['rho3'] + mod5a$coefficients['rho4']
mod5a$urt <- lht(mod5a, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod5a$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod5a$bg <- bgtest(mod5a)$p.value



## Model 6
# change rho to haiti supply price
pri3$rho1 <- pri3$price_sup_ha_log_lag1; pri3$rho2 <- pri3$price_sup_ha_log_lag2; pri3$rho3 <- pri3$price_sup_ha_log_lag3; pri3$rho4 <- pri3$price_sup_ha_log_lag4
mod6 <- lm(price_sup_ha_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             product + date3 - 1,
           data = pri3)
# clustered SEs
mod6$se <- coeftest(mod6, vcov = cluster.vcov(mod6, pri3$product))[,2]
mod6$vcov1 <- cluster.vcov(mod6, pri3$product)
# wild bootstrap clustered SEs
mod6$vcov2 <- cluster.boot(mod6, pri3$product, R = boot, boot_type = "wild")
mod6$se2 <- coeftest(mod6, mod6$vcov2)['coupXautoc', 2]
mod6$p2 <- coeftest(mod6, mod6$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod6$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod6$jsl <- lht(mod6, 'rho1 + rho2 + rho3 + rho4 = 0', vcov. = mod6$vcov)['Pr(>F)'][2,1]
# test for unit root
mod6$coefficients['rho1'] + mod6$coefficients['rho2'] + mod6$coefficients['rho3'] + mod6$coefficients['rho4']
mod6$urt <- lht(mod6, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod6$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod6$bg <- bgtest(mod6)$p.value


## Model 7
# change rho to world supply price
pri3$rho1 <- pri3$price_sup_wo_log_lag1; pri3$rho2 <- pri3$price_sup_wo_log_lag2; pri3$rho3 <- pri3$price_sup_wo_log_lag3; pri3$rho4 <- pri3$price_sup_wo_log_lag4
mod7 <- lm(price_sup_wo_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             product + date3 - 1,
           data = pri3)
# clustered SEs
mod7$se <- coeftest(mod7, vcov = cluster.vcov(mod7, pri3$product))[,2]
mod7$vcov1 <- cluster.vcov(mod7, pri3$product)
# wild bootstrap clustered SEs
mod7$vcov2 <- cluster.boot(mod7, pri3$product, R = boot, boot_type = "wild")
mod7$se2 <- coeftest(mod7, mod7$vcov2)['coupXautoc', 2]
mod7$p2 <- coeftest(mod7, mod7$vcov2)['coupXautoc', 4]
# mean dep var
inc.obs <- complete.cases(pri3[,colnames(mod7$model)])
mean(pri3$coup[inc.obs==T])
# number clusters
length(unique(pri3$product[inc.obs==T]))
# test for joint signif of lags
mod7$jsl <- lht(mod7, 'rho1 + rho2 + rho3 + rho4 = 0', mod7 = mod7$vcov)['Pr(>F)'][2,1]
# test for unit root
mod7$coefficients['rho1'] + mod7$coefficients['rho2'] + mod7$coefficients['rho3'] + mod7$coefficients['rho4']
mod7$urt <- lht(mod7, 'rho1 + rho2 + rho3 + rho4 = 1', vcov. = mod7$vcov)['Pr(>F)'][2,1]
# test for autocorr in error
mod7$bg <- bgtest(mod7)$p.value


## for paper

stargazer(mod1, mod2, mod3, mod4, mod5, mod5a, mod6, mod7, 
          omit = c('date3', 'product', 'Xconf'),
          covariate.labels = c('Coup $\\times$ Autocracy', 'Coup $\\times$ Quake', 
                               'Centrality $\\times$ Autocracy', 'Centrality $\\times$ Quake', 
                               '$\\rho_1$', '$\\rho_2$', '$\\rho_3$', '$\\rho_4$',
                               'World Supply Price', 
                               'Number Firms $\\times$ Autocracy', 'Consumption Share $\\times$ Autocracy'),
          se = list(mod1$se, mod2$se, mod3$se, mod4$se, mod5$se, mod5a$se, mod6$se, mod7$se),
          keep.stat = c('n', 'rsq'), no.space = T, digits = 3, 
          add.lines = list(c('Extra row of SEs on Coup X Autoc', paste0('[',c(round(c(mod1$se2[1], mod2$se2[1], mod3$se2[1], mod4$se2[1], mod5$se2[1], mod5a$se2[1], mod6$se2[1], mod7$se2[1]),
                                                                       digits = 3)),']')),
                           c('Number of Lagged DVs', '', '1', rep('4',6)),
                           c('Month FE', rep('\\checkmark',8)),
                           c('Product FE', rep('\\checkmark',8)),
                           c('Product $\\times$ Conflict Events', rep('',4), '\\checkmark', '\\checkmark',  '', ''),
                           c('Breusch-Godfrey test (p-value)', round(c(mod1$bg, mod2$bg, mod3$bg, mod4$bg,
                                                                       mod5$bg, mod5a$bg, mod6$bg, mod7$bg), 3)),
                           c('Joint Signif. of Lags (p-value)', '', round(c(mod2$jsl[1], mod3$jsl[1], mod4$jsl[1], 
                                                                        mod5$jsl[1], mod5a$jsl[1], 
                                                                        mod6$jsl[1], mod7$jsl[1]), 3)),
                           c('Unit Root (p-value)', '', round(c(mod2$urt[1], mod3$urt[1], mod4$urt[1], 
                                                            mod4$urt[1], mod5$urt[1], mod5a$urt[1], 
                                                            mod6$urt[1], mod7$urt[1]), 4))))




#####
## calculate long-run change
#####

# calculate mean level of coup partic and one sd above mean
m = mean(pri3$coup); sp1 = mean(pri3$coup) + sd(pri3$coup); sn1 = mean(pri3$coup) - sd(pri3$coup)
r = mean(pri3$coup[pri3$product=="Riz"])

# monthly pct increase in prices for one sd above mean vs. mean
im = sp1 * summary(mod5)$coef['coupXautoc',1] - m * summary(mod5)$coef['coupXautoc',1]
i2s = sp1 * summary(mod5)$coef['coupXautoc',1] - sn1 * summary(mod5)$coef['coupXautoc',1]
im = sp1 * summary(mod5)$coef['price_sup_wo_log',1] - m * summary(mod5)$coef['price_sup_wo_log',1]
i2s = sp1 * summary(mod5)$coef['price_sup_wo_log',1] - sn1 * summary(mod5)$coef['price_sup_wo_log',1]

# calculate long run effect of 0 to 1 on coup partic
c1 <- summary(mod5)$coefficients['coupXautoc',1]
clag1 <- summary(mod5)$coefficients['rho1',1]
clag2 <- summary(mod5)$coefficients['rho2',1]
clag3 <- summary(mod5)$coefficients['rho3',1]
clag4 <- summary(mod5)$coefficients['rho4',1]

lr <- c1 / (1-sum(clag1, clag2, clag3, clag4))
lr

# monthly pct increase for rice vs. mean
irs = (r-m) * summary(mod5)$coef['coupXautoc',1]
irs

# calculate long run effect of mean to one sd above mean on coup partic
(sp1-m) * c1 / (1-sum(clag1, clag2, clag3, clag4))

# calculate long run effect of mean to rice on coup partic
(r-m) * c1 / (1-sum(clag1, clag2, clag3, clag4))
exp((r-m) * c1 / (1-sum(clag1, clag2, clag3, clag4)))

prices2 <- read.csv('01_Data/02_Clean/ihsi_0113_clean.csv')
pd1 = mean(prices2$price[prices2$product=="Riz" & prices2$regime=='d1'], na.rm=T)

pd1 * (1 + (r * c1 / (1-sum(clag1, clag2, clag3, clag4))))
rlong = pd1 * exp(r * c1 / (1-sum(clag1, clag2, clag3, clag4)))
rlong

(rlong - pd1)/pd1



#####
## FIGURE 4: Event study
#####

mod1 <- lm(price_log ~ coup*as.factor(date3) + product - 1,    data = pri3)
t <- summary(mod1)$coef

coefs <- data.frame('date' = unique(pri3$date3), 'coup' = seq(-36,108, 1))

coefs <- data.frame('date'=rownames(t[grep('coup:as.factor(date3)*', rownames(t), value = T),]),
                    'coefs'=t[grep('coup:as.factor(date3)*', rownames(t), value = T),'Estimate'],
                    'ses'=t[grep('coup:as.factor(date3)*', rownames(t), value = T),'Std. Error'])

coefs$up90 <- coefs$coefs + 1.64*coefs$ses
coefs$lo90 <- coefs$coefs - 1.64*coefs$ses
coefs$up95 <- coefs$coefs + 1.96*coefs$ses
coefs$lo95 <- coefs$coefs - 1.96*coefs$ses
coefs$date <- gsub("coup:as.factor(date3)", "", coefs$date, fixed=T)
coefs$date2 <- as.Date(coefs$date)
coefs$coup <- ifelse(coefs$date2 >= as.Date('2004-02-01') & coefs$date2 < as.Date('2006-06-01'), 'Autocracy', 'Democracy')

pdf('03_Figures/autoc_event_plot.pdf')
ggplot(data = coefs, aes(x = date2, y = coefs, color = coup)) +
  geom_hline(aes(yintercept=0), lty=2) +
  geom_point() +
  geom_errorbar(aes(ymin = lo90, ymax = up90)) +
  # geom_vline(aes(xintercept=as.Date('2004-02-01')), color = 'red', lty=2) +
  theme_minimal() +
  theme(text = element_text(size=20), 
        axis.text.x = element_text(angle=45, vjust=0.95, hjust=1),
        legend.position = "bottom", legend.title = element_blank(),
        panel.spacing = unit(0.5, "lines")) +
  scale_color_manual(values=c("red", "black")) +
  labs(x = 'Month', 
       y = 'Coefficient on
Month X Coup Participation')
dev.off()


#####
## TABLE A11: Robustness to imposed autocorrelation coefficients
#####

pri3$price_90 <- pri3$price_log - 0.90*pri3$price_log_lag1
pri3$price_925 <- pri3$price_log - 0.925*pri3$price_log_lag1
pri3$price_95 <- pri3$price_log - 0.95*pri3$price_log_lag1
pri3$price_975 <- pri3$price_log - 0.975*pri3$price_log_lag1
pri3$price_100 <- pri3$price_log - 1*pri3$price_log_lag1

mod1 <- lm(price_90 ~ coupXautoc + coupXquake1 + 
             price_sup_wo_log + 
             num_firms_log:autoc + cshare_log:autoc +
             product:noneliteconflict_log +
             product + date3 - 1, 
           data = pri3)
mod1$se <- coeftest(mod1, vcov = cluster.vcov(mod1, pri3$product))[,2]

mod2 <- lm(price_925 ~ coupXautoc + coupXquake1 + 
             price_sup_wo_log + 
             num_firms_log:autoc + cshare_log:autoc +
             product:noneliteconflict_log +
             product + date3 - 1, 
           data = pri3)
mod2$se <- coeftest(mod2, vcov = cluster.vcov(mod2, pri3$product))[,2]

mod3 <- lm(price_95 ~ coupXautoc + coupXquake1 + 
             price_sup_wo_log + 
             num_firms_log:autoc + cshare_log:autoc +
             product:noneliteconflict_log +
             product + date3 - 1, 
           data = pri3)
mod3$se <- coeftest(mod3, vcov = cluster.vcov(mod3, pri3$product))[,2]

mod4 <- lm(price_975 ~ coupXautoc + coupXquake1 + 
             price_sup_wo_log + 
             num_firms_log:autoc + cshare_log:autoc +
             product:noneliteconflict_log +
             product + date3 - 1, 
           data = pri3)
mod4$se <- coeftest(mod4, vcov = cluster.vcov(mod4, pri3$product))[,2]

mod5 <- lm(price_100 ~ coupXautoc + coupXquake1 + 
             price_sup_wo_log + 
             num_firms_log:autoc + cshare_log:autoc +
             product:noneliteconflict_log +
             product + date3 - 1, 
           data = pri3)
mod5$se <- coeftest(mod5, vcov = cluster.vcov(mod5, pri3$product))[,2]

stargazer(mod1, mod2, mod3, mod4, mod5,
          se = list(mod1$se, mod2$se, mod3$se, mod4$se, mod5$se),
          omit = c('date3', 'product'),
          covariate.labels = c('Coup $\\times$ Autocracy', 'Coup $\\times$ Quake', 
                               'World Supply Price',
                               "Number Firms $\\times$ Autocracy", 
                               "Consumption Share $\\times$ Autocracy"),
          column.labels = c('$\\rho = 0.9$', '$\\rho = 0.925$', '$\\rho = 0.95$', 
                            '$\\rho = 0.975$', '$\\rho = 1$'),
          keep.stat = c('n', 'rsq'), no.space = T, digits = 3)



#####
## TABLE A13: Robustness of price results to controls for product characteristics
#####

pri3$rho1 <- pri3$price_log_lag1; pri3$rho2 <- pri3$price_log_lag2; pri3$rho3 <- pri3$price_log_lag3; pri3$rho4 <- pri3$price_log_lag4

mod1 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             num_firms_log:autoc + cshare_log:autoc +
             price_sup_wo_log + 
             product:noneliteconflict_log + 
             bulk_ln_st:autoc + 
             product + date3 - 1, 
           data = pri3)
mod1$se <- coeftest(mod1, vcov = cluster.vcov(mod1, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod1$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

mod2 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             num_firms_log:autoc + cshare_log:autoc +
             price_sup_wo_log + 
             product:noneliteconflict_log + 
              time_con_st:autoc + 
             product + date3 - 1, 
           data = pri3)
mod2$se <- coeftest(mod2, vcov = cluster.vcov(mod2, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod2$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))


mod3 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             num_firms_log:autoc + cshare_log:autoc +
             price_sup_wo_log + 
             product:noneliteconflict_log + 
              pci_value_st:autoc + 
             product + date3 - 1, 
           data = pri3)
mod3$se <- coeftest(mod1, vcov = cluster.vcov(mod3, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod3$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))


mod4 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             num_firms_log:autoc + cshare_log:autoc +
             price_sup_wo_log + 
             product:noneliteconflict_log + 
             ref_lib:autoc +
             product + date3 - 1, 
           data = pri3)
mod4$se <- coeftest(mod4, vcov = cluster.vcov(mod4, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod4$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))



mod5 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             num_firms_log:autoc + cshare_log:autoc +
             price_sup_wo_log + 
             product:noneliteconflict_log + 
               + divis_ln:autoc +
             product + date3 - 1, 
           data = pri3)
mod5$se <- coeftest(mod5, vcov = cluster.vcov(mod5, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod4$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))



stargazer(mod1, mod2, mod3, mod4, mod5,
          omit = c('date3', 'product', 'rho1', 'rho2', 'rho3', 'rho4'),
          covariate.labels = c('Coup $\\times$ Autocracy', 'Coup $\\times$ Quake',
                               "World Supply Price",
                               "Number Firms $\\times$ Autocracy", "Consumption Share $\\times$ Autocracy",
                               "Bulkiness $\\times$ Autocracy", "Time Sensitivity $\\times$ Autocracy",
                                "Complexity $\\times$ Autocracy", "Ref. Price $\\times$ Autocracy",
                               "Divisibility $\\times$ Autocracy")  ,
          se = list(mod1$se, mod2$se, mod3$se, mod4$se, mod5$se),
          keep.stat = c('n'), no.space = T, digits = 3,
          add.lines = list(c('Product $\\times$ Conflict Events', rep('$\\checkmark$', 5)),
                           c('Lagged Dep. Var.', rep(4, 5)),
                           c('Month FE', rep('$\\checkmark$', 5)),
                           c('Product FE', rep('$\\checkmark$', 5)),
                           c("Clusters", rep(18, 5))))





#####
## TABLE A14: Prices of goods imported by coup participators during autocratic 
##   periods using weights based on measures of data quality
#####

temp <- read.csv('01_Data/03_Outputs/product_0911_retention.csv')
temp <- subset(temp, select = c(product, wgt_retain, qty_retain))

pri3 <- merge(pri3, temp, by = 'product', all.x = T)

pri3$rho1 <- pri3$price_log_lag1; pri3$rho2 <- pri3$price_log_lag2; pri3$rho3 <- pri3$price_log_lag3; pri3$rho4 <- pri3$price_log_lag4

mod0 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             num_firms_log:autoc + 
             cshare_log:autoc +
             product:noneliteconflict_log + 
             product + date3 - 1, 
           data = pri3)
mod0$se <- coeftest(mod0, vcov = cluster.vcov(mod0, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod0$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

pri3[,"(weights)"] <- pri3$qty_pct_fam_id

mod1 <- lm(price_log ~ coupXautoc + coupXquake1 + 
           rho1 + rho2 + rho3 + rho4 + 
           price_sup_wo_log + 
           num_firms_log:autoc + 
           cshare_log:autoc +
             product:noneliteconflict_log + 
             product + date3 - 1, 
           data = pri3, weights = qty_pct_fam_id)
mod1$se <- coeftest(mod1, vcov = cluster.vcov(mod1, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod1$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

pri3[,"(weights)"] <- pri3$qty_pct_firm_id

mod2 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             num_firms_log:autoc + 
             cshare_log:autoc +
             product:noneliteconflict_log + 
             product + date3 - 1, 
           data = pri3, weights = qty_pct_firm_id)
mod2$se <- coeftest(mod2, vcov = cluster.vcov(mod2, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod2$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

pri3[,"(weights)"] <- pri3$qty_retain

mod3 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             num_firms_log:autoc + 
             cshare_log:autoc +
             product:noneliteconflict_log + 
             product + date3 - 1, 
           data = pri3, weights = qty_retain)
mod3$se <- coeftest(mod3, vcov = cluster.vcov(mod3, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod3$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

pri3$triple_weight <- pri3$qty_retain*pri3$qty_pct_fam_id*pri3$qty_pct_firm_id
pri3[,"(weights)"] <- pri3$triple_weight

mod4 <- lm(price_log ~ coupXautoc + coupXquake1 + 
             rho1 + rho2 + rho3 + rho4 + 
             price_sup_wo_log + 
             num_firms_log:autoc + 
             cshare_log:autoc +
             product:noneliteconflict_log + 
             product + date3 - 1, 
           data = pri3, weights = triple_weight)
mod4$se <- coeftest(mod4, vcov = cluster.vcov(mod4, pri3$product))[,2]
inc.obs <- complete.cases(pri3[,colnames(mod4$model)])
mean(pri3$coup[inc.obs==T])
length(unique(pri3$product[inc.obs==T]))

## for paper

stargazer(mod0, mod1, mod2, mod3, mod4,
          omit = c('date3', 'product'),
          covariate.labels = c('Coup $\\times$ Autocracy', 'Coup $\\times$ Quake', 
                               '$\\rho_1$', '$\\rho_2$', '$\\rho_3$', '$\\rho_4$',
                               'World Supply Price', 
                               'Number Firms $\\times$ Autocracy', 'Consumption Share $\\times$ Autocracy'),
          se = list(mod0$se, mod1$se, mod2$se, mod3$se, mod4$se),
          keep.stat = c('n', 'rsq'), no.space = T, digits = 3)




#####
## FIGURE A8: Placebo: Coefficient on Coup X Autocracy after arbitrarily moving the window of autocracy
#####

pri4 <- pri3
pri4 <- pri4[order(pri4$product, pri4$date2),]

lags = seq(-48, 48, 1)
coefs = matrix(NA, length(lags), 6)

for (i in 1:length(lags)){
  
  pri4$autoc_lag <- unlist(tapply(pri4$autoc, INDEX = pri4$product, FUN = function(x) shift(x, lags[i])))
  pri4$coupXautoc_lag <- pri4$coup*pri4$autoc_lag
  
  mod <- lm(price_log ~ coupXautoc_lag +
              product + date3,
            data = pri4)
  
  coefs[i,1] <- summary(mod)$coef['coupXautoc_lag',1]
  coefs[i,2] <- coeftest(mod, vcov = cluster.vcov(mod, pri4$product))['coupXautoc_lag',2]
    # coeftest.cluster(pri4, mod, 'product')['coupXautoc_lag',2]
  
}

coefs <- data.frame(coefs, 'lag' = lags*(-1))
coefs[,3] <- coefs[,1] + 1.64*coefs[,2]
coefs[,4] <- coefs[,1] - 1.64*coefs[,2]
coefs[,5] <- coefs[,1] + 1.96*coefs[,2]
coefs[,6] <- coefs[,1] - 1.96*coefs[,2]
coefs <- coefs[order(coefs$lag),]

pdf('03_Figures/autoc_window_nocontrols.pdf')
par(mfrow=c(1,1))
plot(x = NULL, y = NULL, xaxt = "n", xlab = 'Lagged Autocratic Window (Months)', ylab = 'Coefficient',
     #      main = "Placebo Test: Shifting the Window of Autocracy", 
     cex.axis = 1.5, cex.lab = 1.5,
     ylim = c(min(coefs[,1:6], na.rm=T), max(coefs[,1:6],na.rm=T)), xlim = c(min(coefs$lag), max(coefs$lag)))
points(x = coefs$lag, y = coefs[,1], type = "l", col = 'red')
polygon(x = c(coefs$lag, rev(coefs$lag)), y = c(coefs[,3], rev(coefs[,4])), density = NULL,
        col = '#FF003322', border = NA)
abline(h = 0, lty = 2, lwd = 2)
abline(v = 0, lty = 2, col = 'red', lwd = 2)
axis(side = 1, at = seq(-48, 48, 12), seq(-48, 48, 12), cex.axis = 1.5)
abline(v = -26, lty = 2, col = 'black', lwd = 2)
legend('bottomright', lty = c(2, 2), col = c('black', 'red'), lwd = 2,
       legend = c('failed coup', 'real coup'), cex = 1.5)
dev.off()



#####
## FIGURE A9: Robustness: Coefficient on Coup X Autocracy dropping each product
#####


for (i in 1:length(unique(pri3$product))){
  temp <- subset(pri3, pri3$product!=unique(pri3$product)[i])
  temp2 <- lm(price_log ~ coupXautoc + coupXquake1 + price_sup_wo_log + 
                price_log_lag1 + price_log_lag2 + price_log_lag3 + price_log_lag4 + 
                num_firms_log:autoc + cshare_log:autoc +
                product:noneliteconflict_log +
                date3 + product - 1, 
              data = temp)
  temp2$se <- coeftest(temp2, vcov = cluster.vcov(temp2, temp$product))[,2]
  # temp2$se <- coeftest.cluster(temp, temp2, "product")[,2]
  assign(paste0('mod',i), temp2)
}

models <- list(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, mod11, mod12, 
               mod13, mod14, mod15, mod16, mod17, mod18)


pdf('03_Figures/robust_drop.pdf')

# add the data
est <- NA; se <- NA
for (i in 1:length(models)){
  mod <- models[[i]]
  est[i] <- mod$coefficients['coupXautoc']
  se[i] <- mod$se['coupXautoc']
  if (i==21){
    se[i] <- summary(mod)$coefficients[2,2]
  }
}

col = heat.colors(length(unique(pri3$product)))

# create an empty plot for total customization
plot(NULL, ylim = c(-0.02,0.04), xlim = c(0, length(models) + .3), 
     axes = F, xlab = NA, ylab = NA)       

for (i in 1:length(est)) {
  points(y = est[i], x = i, pch = 19, cex = 1, col = col[i])
  lines(y = c(est[i] + 1.96*se[i], est[i] - 1.96*se[i]), x = c(i, i), lwd = 3, col = col[i])
  lines(y = c(est[i] + 1.64*se[i], est[i] - 1.64*se[i]), x = c(i, i), lwd = 5, col = col[i])
}

# add axes and labels
axis(side = 2)      
abline(h = 0, lty = 2, col = "black", xpd=F)
mtext(side = 2, "Coefficient", line = 3)
# mtext(side = 3, "Coefficient on Coup X Autoc dropping each product", line = 1)
legend('bottom', 
       legend = c('Cigarettes', 'Cola', 'Edible oil', 'Kerosene', 'Evap. milk', 'Corn meal',
                  'Medicine', 'Furniture', 'Bread', 'Dry beans', 'Fish', 'Chicken', 'Rice',
                  'Sandals', 'Washing powder', 'Beauty care', 'Sugar', 'Fabric'),
       pch = rep(16, 22),
       col = col, ncol = 3, cex = 1)
box()

dev.off()


